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Abstract 

All attempts of aeroelastic explanations for the torsional instability of suspension bridges have 
been somehow criticised and none of them is unanimously accepted by the scientific community. 
We suggest a new nonlinear model for a suspension bridge and we perform numerical experiments 
with the parameters corresponding to the collapsed Tacoma Narrows Bridge. We show that the 
thresholds of instability are in line with those observed the day of the collapse. Our analysis enables 
us to give a new explanation for the torsional instability, only based on the nonlinear behavior of 
the structure. 
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1 Introduction 

Many suspension bridges manifested aerodynamic instability and uncontrolled oscillations leading to 
collapses, see e.g. mm- These accidents are due to several different causes and the focus of this 
paper is to analyse those due to wide torsional oscillations. Thanks to the videos available on the web 
m, many people have seen the spectacular collapse of the Tacoma Narrows Bridge (TNB) occurred 
in 1940. The torsional oscillations were considered the main cause of the collapse [21 [40]. But the 
appearance of torsional oscillations is not an isolated event occurred only at the TNB: among others, 
we mention the collapse of the Brighton Chain Pier in 1836, the collapse of the Wheeling Suspension 
Bridge in 1854, the collapse of the Matukituki Suspension Footbridge in 1977. We refer to [T8| for a 
detailed description of these collapses. 

These accidents raised some fundamental questions of deep interest for both engineers and math¬ 
ematicians. Due to the vortex shedding, longitudinal oscillations are to be expected in suspension 
bridges, but the reason of the sudden transition from longitudinal to torsional oscillations is less clear. 
So far, no unanimously accepted response to this question has been found. Most attempts of ex¬ 
planations are based on aeroelastic effects such as the frequency of the vortex shedding, parametric 
resonance, and flutter theory. The purpose of the present paper is to show that the origin of the 
torsional instability is purely structural. 

Von Karman, a member of the Board appointed for the Report [2|, was convinced that the torsional 
motion seen on the day of the collapse was due to the vortex shedding that amplified the already 
present torsional oscillations and caused the center span to violent twist until the collapse, see na 
p.31]. But Scanlan [39l p.841] proved that the frequency of the torsional mode had nothing to do with 
the natural frequency of the shed vortices following the von Karman vortex pattern: the calculated 
frequency of a vortex caused by a 68 km/h wind is 1 Hz, whereas the frequency of the torsional 
oscillations measured by Farquharson (an engineer witness of the TNB collapse, the man escaping in 
the video m) was 0.2 Hz , see [TOl p.l20]. The conclusion in m p.l22] is that the vortex trail is a 
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consequence, not a primary cause of the torsional oscillation. Also Green-Unruh [20l § III] believe that 
vortices form independently of the motion and are not responsible for the catastrophic oscillations of 
the TNB. The vortex theory was later rediscussed by Larsen [27l p.247], who stated that vortices may 
only cause limited torsional oscillations, but cannot be held responsible for divergent large-amplitude 
torsional oscillations. Recently, McKenna [33] noticed that the behavior described in Larsen’s paper 
was never observed at the Tacoma Bridge and also Green-Unruh [20] believe that the Larsen model 
does not adequately explain data or simulations at around 23m/s. 

Bleich m suggested a possible connection between the instability in suspension bridges and the 
flutter speed of aircraft wings; but Billah-Scanlan uni p.l22] believe that it is a great mistake to 
relate these two phenomena. Billah-Scanlan also claim that their own work proves that the failure of 
the TNB was in fact related to an aerodynamically induced condition of self-excitation in a torsional 
degree of freedom; but Larsen [271 P-244] believes that the work in [10] fails to connect the vortex 
pattern to the switch of damping from positive to negative. Moreover, McKenna [33] states that m is 
a perfeetly good explanation of something that was never observed, namely small torsional oseillations, 
and no explanation of what really oeeurred, namely large vertieal oseillations followed by torsional 
oseillations. 

The parametric resonance method was adapted to a TNB model by Pittel-Yakubovich [35l [36] , see 
also @31 Chapter VI] for the English translation and a more general setting. The conclusion on 
[43]. p.457] claims that the most dangerous phenomenon for the stability of suspension bridges is a 
eombination of parametrie resonanee. But Scanlan [39l p.841] comments these attempts by writing 
that Others have added to the eonfusion. A reeent mathematies text for example, seeking an 
applieation for a developed theory of parametrie resonanee, attempts to explain the Taeoma Narrows 
failure through this phenomenon. We refer to @21 ED El for connections between (aerodynamic) 
parametric resonance, flutter theory, and self-oscillations, also applied to suspension bridges. 

To conclude this quick survey of attempts for aeroelastic explanations, we mention that Scanlan 
[38[ p.209] writes that ...the original Taeoma Narrows Bridge withstood random buffeting for some 
hours with relatively little harm until some fortuitous eondition ^Lroke’^ the bridge aetion over into its 
low antisymmetrieal torsion flutter mode] the words fortuitous eondition tell us that no satisfactory 
aeroelastic explanation is available. Due to all these controversial discussions, McKenna [32l § 2.3] 
writes that there is no eonsensus on what eaused the sudden ehange to torsional motion, whereas Scott 
@01 writes that opinion on the exaet eause of the Taeoma Narrows Bridge eollapse is even today not 
unanimously shared. Summarizing, all the attempts to find a purely aeroelastic explanation of the 
TNB collapse fail either because the quantitative parameters do not fit the theoretical explanations 
or because the experiments in wind tunnels do not confirm the underlying theory. 

Nowadays the attention has turned to the nonlinear behavior of structures [3l[26]. In a recent paper [4] 
we gave an explanation in terms of a struetural instability: we considered an isolated nonlinear bridge 
model and we were able to show that, if the longitudinal oscillations are sufficiently small, then they are 
stable whereas if they are larger they can instantaneously switch to destructive torsional oscillations. 
The main tool used there are transfer maps (Poincare maps), which highlight an instability when the 
characteristic multipliers exit the complex unit circle. The same phenomenon was later emphasized 
for different models using the instability for the Hill equation, see IDEIEI, which can also be explained 
by the Floquet theory. All these results were obtained by considering isolated systems, that is, by 
neglecting both the aerodynamic forces and the dissipation. The idea to consider an isolated system 
was already suggested by Irvine [23l p.l76] for a suspension bridge model similar to the one considered 
in [1]: Irvine ignores damping of both structural and aerodynamic origin, his purpose being to simplify 
as much as possible the model by maintaining its essence, that is, the conceptual design of bridges. 
And since our purpose is precisely to highlight the role of the structure on the stability, we follow this 
suggestion: in the conclusions of the present paper we mention how the aerodynamic effects combine 
with the structural effects. 

In this paper we improve the results in [4] by introducing a more precise model which also takes 
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into account the nonlinear restoring action of the cables+hangers system on the deck; the model is 
described in detail in Section [2] and is the nonlinear version of a linear model that we introduced in 
[5]. In Section]^ we define the nonlinear longitudinal modes, that is, the periodic purely longitudinal 
motions of the deck which may be unstable and create torsional motions. In Sections |4.1| and |4.2| 
we describe the theoretical framework that we use to study the stability of the longitudinal modes 
and we prove that they are stable for small energies. Finally, we validate our theory by performing 
numerical simulations of the model with the parameters of the collapsed TNB. We study the existence 
and the behavior of approximate longitudinal modes and we compute their instability thresholds. 
Our numerical results confirm that the longitudinal oscillations with 8 or 9 nodes and amplitudes of 
about 4m, namely the oscillations observed the day of the collapse, are prone to generate torsional 
oscillations, see Section 4^ The main conclusions and our explanation of the TNB collapse are 
collected in Section [H 


2 The mathematical model 

2.1 Description of the structure 

Throughout this paper the variable x indicates the position along the deck and we denote the space 
and time derivatives as follows: 

/X f dw , . , dw , . I ^ It 

w — wix) ^ w =, w — wit) w — , w — w[xA) ^ < 

^ ^ dx' ^ ^ dt ' ^ ^ \ ^ 

{ ~ dt 

and similarly for higher order derivatives. The physical constants are listed in Section [4.4[ 

In a suspension bridge four towers sustain two cables that, in turn, sustain the hangers. At their 
lower endpoint the hangers are linked to the deck and sustain it from above. The hangers are hooked 
to the cables and the deck is hooked to the hangers. In this section we quickly revisit the bridge 
model that we recently introduced in [5]: the displacements involved were assumed to be small, which 
justifies the asymptotic expansion of the energies. Here we refine this model by considering more 
terms: this raises several nonlinearities in the Euler-Lagrange equations. Moreover, we neglect the 
flexibility of the hangers: this will be justified below. As usual in engineering literature, we assume 
that positive displacements of the deck and the cables are oriented downwards and the origin is at the 
level of the deck at rest. 

We assume that the deck has length L and width 2i with 2£ <C T; we model it as a degenerate plate 
consisting of a beam representing the midline and cross sections which can rotate around the beam. 
The beam contains the barycenters of the cross sections and we denote its position with y — y{x^t). 
The angle of rotation of the cross sections with respect to the horizontal position is denoted by 
9 — 9{x^ t). Then the positions of the free edges of the deck are given by ^ ± £ sin 9 ^ y ±19 since we 
aim to study what happens in a small torsional regime. We emphasise that our results do not aim to 
describe the behavior of the bridge when the torsional angle becomes large; instead, we explain how 
a small torsional angle can suddenly increase, that is, our purpose is to describe the mechanism that 
triggers torsional oscillations. 

Since the spacing between hangers is small relative to the span, the hangers can be considered as a 
continuous membrane connecting the cables and the deck. We denote by —s{x) the position of the 
cables at rest, —sq < 0 being the level of the left and right endpoints of the cables (sq is the height of 
the towers); L is the distance between the towers. If a beam of length L and linear density of mass 
M is hanged to a cable whose linear density of mass is m, and if g denotes the gravitational constant, 
then the cable at rest is subject to a downwards vertical force density given by (M + my^l + s'{x)‘^) g. 
Each free side of the deck is connected with hangers to a cable so that each cable sustains the weight 
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of half deck; then, we obtain the following boundary value problem for s(x): 



It is not difficult to show that Q admits a unique solution which, moreover, is symmetric with respect 
to X = i^/2, see [SJ Proposition 1] for the details. In order to simplify notations we set 

^(x) := a/ 1 + s'(x)2 , 

SO that, ^(x) represents the local length of the cable at rest. When the deck is mounted, the hangers 
are in tension and reach the length s(x), where s is the solution of 0 ; if no additional load acts on the 
system, the equilibrium position of the deck is horizontal and the position of the cables at equilibrium 
(deck mounted) is —s(x) for x E (0,L). 

The two sustaining cables are labeled with i = 1,2. If pi = pi{x^t) denotes the displacement of the 
cable i with respect to the rest position, then pi{x^t) — s(x) denotes the actual position of the cable. 

The flexibility of the hangers has a significant effect on the frequencies of the higher modes when the 
deck is very stiff, see m- However, as far as lower modes and weakly stiffened bridges are involved, 
the elastic deformation of the hangers may be neglected, see again the results by Luco-Turmo m- 
Since we aim to model a very flexible deck as the TNB, we assume that the hangers are rigid, their 
action merely being to connect the deck to the cables. This assumption is justified by the fact that the 
action of the cables is the main cause of the nonlinearity of the restoring force, see Bartoli-Spinelli isi 
p.180]. Moreover, the nonlinear contribution of the hangers is mainly due to their slackening but, as 
reported in [2], the slackening of the hangers occurred at the TNB only after that the large torsional 
oscillations appeared. Our purpose is to understand how negligible torsional oscillations of the deck 
suddenly become dangerous ones, that is, to describe what happens before the slackening starts. This 
is a further reason to consider rigid hangers. If we assume that the hangers have a fixed length, then 
Pi — y -\- W and P2 = y — ^0: for our convenience, in the next subsection we maintain the notation 
with the pfs while we replace them in the Euler-Lagrange equations 0-0. 

2.2 Kinetic and potential energies of the structure 

The kinetic energy of the deck is the sum of the kinetic energy of the barycenter of the cross section 
and of the kinetic energy of the torsional angle. Then the total kinetic energy of the bridge is 

Ekin = Y Jo iP2)tMx) dx . 

The gravitational energies of the deck and the cables are, respectively, given by 

pL pL 

—Mg / ydx and —mg / [pi-\-P 2 )i{x) dx . 

Jo Jo 

The elastic energy of the deck is composed by the bending energy of the beam and the torsional 
energy: 



The tension of each cable consists of two parts. The first part is the tension at rest H{x) = Hq^{x). 
The amount of energy needed to deform each cable at rest under the tension H(x) in the infinitesimal 
interval [x^x + dx] from the original position —s{x) to —s{x) -\- Pi{x^t) {i = 1,2) is the variation of 
length times the tension, that is 

Ec{x) dx = Hq^{x) (a/1 -I- [s'(a;) - {pi)^{x,t)]'^ - y/l + s'{xy^ dx. 
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Hence, using the asymptotic expansion 


- yr+w= 


s'iP: 




^ ' 2e 2^5 

and neglecting the higher order terms, the energy necessary to deform the whole cable is 

ErXPi) = J Ec{x) dxm Ho J (^-s'{pi)x + 1 dx . 


( 2 ) 


( 3 ) 


2^2 2e 

The second part of the tension of each cable is the additional tension 

AE^, , 

due to the increment of length r(p^) of the cable (i = 1,2). This increment requires the energy 

A F' 

EtciPi) = ^ . 

In view of the asymptotic expansion Q, we can approximate the increment r(p^) of the length of the 
cable (due to its displacement pi) by 


( 4 ) 


r(Pi) J yi + [s' - {Pi)x? dx- Lc^ J ^ 


+ 


2e 


A third order expansion then gives 




-s'{P i)x , {Pi) 

{Pi)x 


dx . 


/o C 

Therefore, we approximate @ with 


Etc{Pi) — 


s' 

Jo 


s'ijH 


AE 


rT)(r 




iJx 


( 5 ) 


2.3 The Euler-Lagrange equations 

Some integration by parts show that the variation of the energy Ec{pi) in Q is given by 

{Pi 


dEc{pi)z = Hq 


s — 


3s'{pX‘^ 


tJx 


'0 V e 2e 


2 ; dx. 


For the variation of the energy Etc{pi) in Q, we neglect terms of order larger than 2, and by integrating 
by parts we obtain 


dEfci^Pi) z — 


AE 


rM^ r 

Jo e JJo 


s" . AE 
^ zdx + —— 
e Lc 


ps"pA f 

Jo e J Jo 


L 


(pi) 


2 ; dx . 


By recalling that pi = y + W and p 2 — y — l9^ the first variation of the sum of all the energies involved 
yields the following Euler-Lagrange equations: 


{M+2m^)yu = -Ely x: 

2AE 


+ Ho ^ + 3 


s"y 


Lc Vo 
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This is the system describing the dynamics of a suspension bridge that we will consider 
throughout the paper. Since the degenerate plate is hinged between the two towers and the cross 
sections between the towers cannot rotate, the boundary conditions to be associated to read 

y{0,t) = y{L,t) = yxx{0,t) = ya;x{L,t) = 9{0,t) = e{L,t) = 0 Vt > 0 . (8) 


For our investigation of (©-0 (both theoretical and numerical) we first set 9 = 0 in a to obtain 

(9) 


'.,2 

X 


‘2iU u 

{M-\-2m^)ytt-\-EIyxxxx = Hq ( ^ 


AE 

\f^yl] 

s" 

2AE 

cce" 
_1 
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[Jo e\ 
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Lc 

[Jo e\ 
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and we find periodic solutions y of different amplitude and number of zeros. Then, for each solution 
y of 0 we consider 0 as a linear equation in 9 with periodic coefficients and we study the stability 
of the trivial solution of 0. 

Our characterization of torsional stability consists in taking a periodic solution y of Q and plugging 
it into 0; we say that y is torsionally stable if the trivial solution of the so obtained equation (see 
( [T^ below) is stable. This is basically the characterization suggested by Dickey [16] for a couple of 
oscillators related to a nonlinear string: it is also called linear stability in literature [191 Definition 2.3]. 
An apparently stronger form of stability, without linearization and substitution of ?/, is the well-known 
Lyapunov characterization. For a simpler system of ODEs Ghisi-Gobbino m were able to prove that 
these two characterizations are equivalent. 


3 Nonlinear longitudinal modes 

The milestone paper by Rabinowitz m started the systematic study of the existence of periodic 
solutions for the nonlinear string equation 


~ "^xx + /(^5 u) = 0 (x, t) G (0, tt) X , ?i(0, t) — it(7r, t) = 0 t G . (10) 

Prior to m, only perturbation techniques were available and the nonlinearity was assumed to be small 
in a suitable sense. Rabinowitz proved that (10) admits periodic solutions. In fact, there exist infinitely 
many periodic solutions of (10), having “almost” any possible period T > 0: quite sophisticated tools 
are necessary to obtain a full description of the set of periodic solutions, see e.g. nzi and references 
therein. Periodic solutions were also found for the related nonlinear beam equation 


T "^xxxx T f (^5 '^) — 0 (x, t) G (0, tt) X , 

7i(0, t) = u{7T, t) = Uxx{^, t) = Uxx{ 7T, t) = 0 t G M+ , 


( 11 ) 


see 


[23130]. Also nonlocal (and nonlinear) variants of ( [Tol ) have been considered. In 1876 Kirchhoff 
introduced the following nonlinear vibrating string equation 


utt- (a + 


(^a + b J Uxx — 0 (x, t) G (0, tt) x , i^(0, t) — u{'k^ t) = 0 


t G 




( 12 ) 


where a, 6 > 0. Note that (10) with /(x, u) — au + bv? is quite similar to (12) and, in some respect, (12) 
appears more complicated since first order derivatives are involved. Nevertheless, the modal analysis 
of (12) is much simpler since it behaves as a “hidden” linear problem. Exploiting this peculiarity. 


Cazenave-Weissler m proved that ( [l^ has infinitely many periodic solutions. 

If we take 0(x,O) = ^^(^^O) = 0 as initial data in 0-0-0, then 9{x^t) = 0, whereas y solves the 
nonlinear nonlocal equation 0. A careful look at ^ shows that its structure, its nonlinearities, and 
its nonlocal terms are of the same kind as those appearing in and ( [T^ . It is therefore 

reasonable to conjecture that 0 admits infinitely many periodic-in-time solutions. 
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In order to find such solutions, or at least a reliable approximation of them, we proceed numerically. 
The first step is to look for small periodic solutions of © close to those of the linearized problem 


(M + 2m^)ytt + Ely XXXX — 2Ho ( p 


2AE 


J 


^ s"y 




* 


(13) 


Let us analyze (13) and see what has to be expected from the numerical procedure. Periodic solutions 
of (13) may be obtained by separating variables. Consider the linear operator C defined by 

.f\ ' 

^ 


Cu EIu”” -2 Hq[^] + 


u 
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L^o 


^ Vn€C^[0,L]. 


The definition of £ may be continuously extended to the Sobolev space fl Hq{0,L): the operator 
£ maps this space into its dual space H' and is implicitly defined by 


{Cu, v) ^ El [ 
Jo 


u'E 


u"v" + 2Ho [ + 

Jo e 


2AE 

■ p s''u 


■ p s''v 

Lc 

Jo 'W. 


Jo 'W. 


^ yu,veH^n (0, L) (14) 


where (•, •) denotes the duality pairing between n Hq{{),L) and H'. The bilinear form in (14) is 
symmetric and coercive. Therefore, the operator C is self-adjoint and the eigenvalue problem 


Cu — \{M + 2m^)u 


(15) 


admits an unbounded sequence {Xk} of positive eigenvalues whose corresponding eigenfunctions form 
a complete system in H (0, L). 

If the cables were in the limit horizontal position, that is s' — s" — {) and ^ = 1, then the eigenfunc¬ 
tions of (jlSj would be sin(^x) for all integer A: > 1 with corresponding eigenvalues given by 


(M -h 2m)L^ 


Xk = 


As we could verify numerically, when the shape of the cable satisfies the eigenfunctions of ( [T^ 
remain close to sin(^x) with A: —1 zeros in (0,L). 

Denote by (px an eigenfunction of (15) relative to the eigenvalue A, then for any A, S G M the function 


u{x^ t) = pj\{x) ^Asin (V^t) -h B cos {VXt)^ 


(16) 


is a periodic solution of (13). 

We numerically find approximate periodic solutions for ^ in two steps. First, we transform the PDF 
into a system of ODEs by discretising in space via Fourier polynomials and we build a Mathematica 
program that solves the corresponding initial value problems. Then we look for points in the (finite 
dimensional) phase space that belong to a periodic orbit. 

More precisely, given the boundary conditions ([^, we look for an approximate solution of ^ that 
can be represented by a trigonometric polynomial 


y{x,t) = ^yfc(t) sin 


k=l 


/ klTX 


(17) 


To achieve this goal, we plug 0 into © and we project it onto the span of {sin (ttx/L) ,..., sin {nnx/L)}: 
in particular, we choose initial conditions 


k=i 


knx 


y{x, 0) = yok sin I —— ) , yt{x, 0) = yik sin 


/ knx 


k=l 


V L 
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for some Yq = {yok}k=i,...,n ^ sind Yi = {yik}k=i,...,n ^ The dimension n is chosen in such 
a way that the result of the numerical integration does not change significantly when n is increased. 
For the computation of the first 6 periodic solutions we take n — 10, while for the computation of the 
next 4 periodic solutions we take n = 16. 

Let Y{t) = {yk{t)}k=i,...,n ^ we obtain a system of n nonlinear ODEs in the form 

Y{t) = G{Y{t)) {G:R^ ^ M^) (18) 


with initial conditions T(0) = Yq and T(0) = Yi. We use Newton’s algorithm to find initial vectors 
lo and Yi that lead to periodic solutions of ( [T8| ). For any T > 0 we define the transfer map 

$r : ^ , $r(Yo, Yi) = (Y(r), Y(r)) , 


where Y{t) is the solution of (18) with initial conditions (Yq, Yi), and we look for fixed points of 
by iterating the transformation 


N{Yo,Yi) (yo,n) - [J^TiYo,Yi) - I]-\^TiYo,Yi) - iYo,Yi)) 

where J^t is the Jacobian matrix of and / is the identity 2n x 2n matrix. As initial point for the 
Newton’s algorithm we choose Yq = ack and Yi = 0, where a is some positive number, G {1,..., n} 
and {ei,...,e^} is the canonical basis in We observe that, if a is sufficiently small, there exists 
T > 0 such that $r(Yo, Yi) (Yq, Yi): for such T, Newton’s algorithm converges rapidly. In Table 
we quote the period T of each mode obtained for cr ^ 0. 


Mode 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Period (s) 

10.95 

7.67 

5.42 

3.75 

2.9 

2.41 

2.02 

1.72 

1.5 

1.32 


Table 1: Period of the longitudinal mode on each branch for energies close to 0. 


Figure [^represents the graph of the solution obtained by this algorithm fo r k = 8 at time t = 0 


as 


expected, it resembles a multiple of sin(^x) (with L — 853.44m, see Section 
has small L^-norm, we expect it to be an approximate periodic solution of t 


4 . 4 F Since this function 
re linear equation (13). 



Figure 1: y{x^0) for the solution on the 8th branch at T = 1.72 (small energy). 


This small periodic solution of the nonlinear system (18) is an example of what we expect to observe 
when a still bridge is slightly perturbed by some external forcing, typically the wind. If more energy 
is inserted into the structure, the oscillation becomes wider. Since the insertion of energy is gradual, 
it seems natural to look for other solutions, of larger energy and L^-norm, by continuity. To do 
so, we slightly increase T and use again Newton’s algorithm: it turns out that, on each branch, the 
period is an increasing function of the energy and, if the increment in T is sufficiently small, Newton’s 
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algorithm keeps converging rapidly, and the resulting solution of (18) has larger energy and L^-norm 


By increasing T we build the branch of the fc-th nonlinear modes of (18). This procedure cannot be 


extended indefinitely; at some point the convergence of Newton’s algorithm starts deteriorating and 
smaller increments in T must be taken in order to keep following the branch of solutions. We stop the 
procedure when the acceptable increments become too small to be feasible. In this way we obtain all 
the branches of the k-th nonlinear longitudinal modes for A: = 1,..., 10. 


Definition 1. We call a periodic solution of (18) obtained with this algorithm a k-th. nonlinear 
longitudinal mode of period T. 


We observed that the nonlinear longitudinal modes on the odd branches are symmetric with respect 
to the center of the bridge, see e.g. Figure where the solution on the 7th branch at T = 2.18 is 
displayed. On the other hand, the modes on even branches have no symmetries, see e.g. Figure 
where the solution on the 8th branch at T = 1.86 is displayed. 



Figure 2: for the solution on the 7th branch at T = 2.18 (critical threshold of instability). 



Figure 3: ^(x,0) for the solution on the 8th branch at T = 1.86 (critical threshold of instability). 


We remark that the k-th Fourier component of a k-th nonlinear longitudinal mode is much larger 
than the other ones, see e.g. Figure where the evolution in time of the first 12 Fourier components 
for the solution obtained with A: = 8 with T = 1.86 is displayed. 

In all our experiments we also observed that the first Fourier component of the periodic solution 
(namely the function yi{t) in ([T7|)) is negative for all times, see e.g. the first plot in Figure]^ This 
means that the deck is bent upwards during periodic motions. In fact, our results show that, when 
the periodic motion starts, the cables undergo an additional tension due to the nonlocal term and 
therefore they pull the deck upwards. 


9 











0.D010 

0.D0D5 


-0.D0D5 

-0.D010 

0.D04 

0 .D02 


- 0 .D 02 

-0.D04 



O.DOG 

-D.D04 

0.D02 


-0.D02 

-0.D04 

-O.DOG 





Figure 4: Time evolution of the Fourier components 1 to 12 of the solution on the 8th branch for 
T = 1.86 (critical threshold of instability). 
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4 Torsional stability of longitudinal modes 


4.1 Theoretical framework 


The purpose of this section is to insert the stability analysis of suspension bridges within a suitable 
theoretical framework. We explain why, for small longitudinal oscillations, one expects the system 
®-0 to remain torsionally stable, and why, for large longitudinal oscillations, the system becomes 
unstable. Most of the tools introduced in this section are applications to our context of classical tools 
from the theory of ODE’s that can be found e.g. in [31131. 

We consider a (time) periodic solution y of ([^ and we denote by Ty its period. Then we view the 
equation 0 for 9 as an independent equation, that is, a linear hyperbolic equation of the form 


a 


{x)9u = (ci{x^t)9x) + ^{x) f b{x^t)9{x^t)dx +c(x,t) f ^{x)9{x^t)dx (19) 


for suitable coefficients a, 6 , c, cr,/3, 7 ; note that the coefficients a, 6 , c depend on the periodic solution 
y = y{x,t) of Q. 

Because of energy conservation, if the energy of the system is at some time concentrated on low 
Fourier coefficients of the solution, e.g. if the initial condition is smooth, then it remains so for all 
times. Hence, it seems sensible to simplify the problem by studying its projection onto the space 
spanned by {sin(Jx),..., sin(^x)}, for some integer v >2. This integer v needs not (so che a volte 
si usa, ma non e corretto. ’’does not need” e meglio) to be the same as the integer n in Section]^ 
in order to state and prove a sufficient condition for stability we simply take v — 2 whereas for our 
numerical experiments we take v — n. 

For each A: = 1, ...,z/ we multiply (19) by sin(^x) and we integrate over (0,L): we obtain z/ ODEs 
in the form 


hit) + ^Xjk{t)0j{t) = 0 , {k = 1,...,0 


( 20 ) 


for some coefficients Xij depending only on t. We take the solution of (20), that is, 

0 (x, t) — 9j {t) sin (^x 


j=i 


( 21 ) 


as an approximate solution of (19). Since the time-dependent coefficients in (19) are computed in 
terms of the periodic solution ^ of ([^, we have that 

the functions t Xjk{t) Ty - periodic Vj, k = 1,..., z/ . (22) 

Moreover, these coefficients are symmetric: Xij{T) = Xji{T) for all i, j. We define the symmetric z/ x z/ 
matrix 


2(t) = (: 


Xjk 


which is Ty - periodic by (22). Then we rewrite 


in the vectorial form 


W + H(t)W = 0, 


(23) 


where W{t) — {9i{t)^ ...^9y{t)) is the finite dimensional approximate vector containing the torsional 
components of the motion. Following 113 p.l55, vol.l] we give the 


Definition 2. We say that the system (23) is strongly stable if all its solutions are hounded in 
and this property is preserved under small perturbations ofS{t). 
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In turn, since the system (23) has coefficients which depend on the particular nonlinear longitudinal 


mode used to build the symmetric matrix S(t), we may characterize the stability of the modes as 
follows. 


Definition 3. We say that a nonlinear longitudinal mode of 


(see Definition^ is torsionally 


stable if the system (23) is strongly stable; otherwise, we say that it is torsionally unstable. 


Clearly, the boundedness of all solutions of the linear system (23) is equivalent to the stability of 


the trivial solution W{t) = 0 and, therefore, to the torsional stability of the nonlinear longitudinal 


mode. Our purpose is precisely to study the stability of (23). In Section 4.2 we state and prove a 
sufficient condition for the stability in the case v — 2 and we discuss how, for general systems such 


as (23), instability may arise. In Section 4.3 we compute numerically the thresholds of instability in 


dependence of the behavior of the nonlinear longitudinal mode which is used to built the coefficient 


matrix S in (23). 


4.2 The case v = 2 


A fairly useful tool to prove the stability for linear systems such as (23) is the following criterion, see 
Theorem II p.270, vol.l]: 


Lemma 1. Assume that a n x n matrix P(t) is eontinuous and T-periodie for some T > 0 and write 
P(t) = D(t) + R(t) where D(t) is the diagonal matrix eontaining the same elements as P(t) and 
R(t) = P(t) — D(t) is the ^^residual matrix’^ eontaining all the eoeffieients off the diagonal o/P(t). 
Let ri and r 2 be two eontinuous T-periodie funetions sueh that 

ri(t)I<R(t) <r 2 (t)I VtG [0,r]. 

Let Sk{t) (k — 1, be the elements of the diagonal matrix D(t). Assume that for all s G [0,1] the 

sealar Hill equations 


'^k + 


+ sri{t) + (1 - 5 


= 0 


{k = 


(24) 


are strongly stable. Then the system Z + P(t)Z = 0 is strongly stable (here, Z — Z{t) ^ 

In this section we study in detail the case v — 2, since this case is simpler to handle and allows to 
prove a sufficient condition for the stability. Moreover, according to the detailed analysis of the TNB 
collapse by Smith-Vincent HD p.21], the only torsional mode whieh developed under wind aetion on 
the bridge or on the model is that with a single node at the eenter of the main span, see also [18] for 
further evidence showing the importance of the second torsional component. 


If z/ = 2 then (23) simply becomes 



Xn{t) 0 

0 X22it) 



0 Xi2{t) 

Xuit) 0 




(25) 


where we have emphasized the diagonal matrix D(t) and the residual matrix R(t). We point out that 
a simple and elegant form such as ( [2^ is made possible by the fact that n — 2: the eigenvalues of the 
residual (symmetric) matrix R(t) are ±xi 2 (t) which means that 


-|xi 2 (t)|i<R(t)<|xi 2 (t)|i vtG[0,r,]. 


By (26), in the situation under study the n — 2 families of Hill equations in (24) read 

h + (xii(i) + OiXi2{t)fi = 0 
^2 + (x22{t) + OiX\-2{t)]^‘2 — 0 


(-1 < a < 1 ). 


(26) 


(27) 


Let us derive (27) in the limit case where y = 0. 
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Lemma 2. There exist 71,72 > 0 sueh that, if y = 0 , then the funetion 

+ 02 {t) sin 

satisfies the projeetion of 0 onto span{sin( Ja:), sin(^a:)} with Oi and 62 solving 

01+7101 = 0 , 02 + 7202 = 0 . ( 29 ) 


9{x, t) = 0 i{t) sin 



Proof. If y = 0, then ([^ reads 


M 


— + 2mn eOtt = + 2^1 7 - ^ 


0 . 


AE 


/o 


e. 


s" 

e 


(30) 


and all the coefficients are now independent of t. The functions involved in this equation are either 
symmetric or skew-symmetric with respect to x = ^: 

s"{cc),^{cc),smGx),cos(^x) are symmetric, s'{cc),sm0rx),cos{j-x) are skew-symmetric. (31) 


Several terms cance 


Take 0 in the form (28) and plug it into (p0|; then multiply (30) by sin(TaO and integrate over (0, L). 


due to (31): if an odd number of skew-symmetric functions is involved then the 


integral over (0, L) vanishes. We end up with the first equation in 


71 = 


1 (GKtt^ 27r‘^fHo 7 cos'^if^) . 

' ' ' ^ dx + 


A\ 2 L 


+ 


L2 


77 


■ s" 

Jo 


with 


s"{x) sin( Jx) 


n 2N 


7a;) 


dx 


> 0 


and 


equation in (29) wit 


A = 


7 J + 2m^{x)^ sin^ ' 


Then multiply (30) by sin(^a:) and integrate over (0, L). By using again (31), we find the second 


72 = 


1 (2GK7r^ Sir^fHo p cos^(^x) 


B 




+ 


L2 


/ 




-dx > 0, B — 




M 


27T 


The proof is complete. 


+ 2 m^{x) j sin^ f -j-x ) dx . 


□ 


Needless to say, the equations (29) admit periodic (bounded!) solutions and 02 for any initial 
data. We denote by N the set of nonnegative integers and we assume that 


TT 


for A: = 1 , 2 


(32) 


where Tq = lim||^||^^o is the period of the corresponding periodic solution of the linear problem. 
Clearly, (32) has probability 1 to occur among all possible choices of {G, K, M,To, L,i) G M^. 

Let Ty be as in (22); if (32) holds, we know that for fc = 1, 2 there exists such that 




2.^2 


<lk< 


{uk + l)^7r' 

rp 2 


( 33 ) 


If ^ ^ 0 and y,yx^ Uxx 0 uniformly with respect to t, then 

Xkk{t) Ik {k — 1, 2), X 12 J) -P 0 


uniformly 
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and the system (27) converges to (29), see Lemma Therefore, by (33) and by continuity, we infer 
that there exists 6 : > 0 such that if || 2 /||oo < ^ then 


7^2 

y 


{nk + l)^7r2 

7^2 

y 


< Xkkit) - 1x12(i)I < Xkkit) + 1x12(i)I < 


for A; = 1,2 and for all t > 0. Then a stability criterion for the Hill equation due to Zhukovskii im, 
see also [l3l Test 1 p.697, voL2], states that the two equations in (27) are strongly stable. In turn. 
Lemmaimplies that the system ( [2^ is strongly stable. We have so proved the following statement. 


Theorem 1. Assume (32) and let y he a nonlinear longitudinal mode of ( [T^ . There exists 6 : > 0 
sueh that i/ H^/Hoo < then the projeetion (23) of 0 onto span{sin( Jx),sin(^x)} is strongly stable; 
therefore, y is torsionally stable. 


Clearly, a solution y of (18) with small L^-norm also has small energy. Theorem says nothing 
about the stability of (23) when the periodic solution y of 0 has large amplitude (or energy). We 
now address this question. 


Instead of (23) we consider a first order linear differential system such as 


Z^A{t)Z, 


(34) 


where Z G for some z/ > 2 and A is a T-periodic 2z/ x 2z/ matrix for some T > 0; the system (23) 
may be written in the form (34) by setting Z = (VL, W). Define a transfer map Z{t) i-G Z{t -\- T) in 


terms of a suitable transition matrix M, that is, Z{t + T) — M Z{t) with Z being a solution of 


It is clear that if all the (complex) eigenvalues of M have modulus 1, then the solutions of (34) are 


globally bounded and therefore the trivial solution Z = 0 of (34) is stable. The numerical procedure 


described in Section computes these eigenvalues. 

Ortega m considers a family of nonlinear Hill equations such as 

^ + a{t)z + f{z) = 0 with a{t -\-T) — a{t) ; 


(35) 


a typical example is f{z) = z^ but many other nonlinearities / yield the same behavior. Taking 
Z := (z, i), this equation may be reduced to a first order nonlinear 2 x 2 system of the kind of ( [M] ): 

Z = Gl(t)Z + F(Z). (36) 

Under suitable conditions on the nonlinearity F, Ortega proves that if the trivial solution of ([34|) is 


stable, then also the trivial solution of (36) is stable. This sufficient condition shows that a linearization 


process is legitimate, if one aims to study the stability of (36). 


In their paper, Cazenave-Weissler m analyzed the stability of simple modes of by using tech¬ 
niques from the Floquet theory, combined with Poincare maps and the Hill equations. They were able 


to prove that, if the energy is sufficiently large, then simple modes of ( 12 ) become unstable: they take 
advantage of the fact that ( 12 ) is “almost linear”. 


For a different nonlinear and local system of ODE’s it was proved in [9] that the same fact occurs, 
large energies yield instability of modes. The analysis in [9] is fairly precise and is made possible by 
the fact that the linearized dynamical system leads to some Mathieu equations (a particular case of 
the Hill equations) for which more precise theoretical tools are available. 

Therefore, both for nonlinear nonlocal problems m and for nonlinear local problems [9], the insta¬ 
bility of the nonlinear longitudinal modes of (18) with large energy has to be expected. The system 


0-0 is certainly much more complicated than the systems considered there, and a rigorous proof of 
the torsional instability at high energies seems out of reach, but we believe that it is governed by the 
same phenomena. In the next section we will see that the numerical experiments support our belief. 
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4.3 Thresholds of torsional instability 


Take a k-th nonlinear longitudinal mode y of (18) having period Ty (see Definition]^ and use it 
to build the coefficient matrix S(t) in (23). Then transform (23) into (34) by setting Z = (W,W). 
We choose 2z/ different initial conditions Z{0) = where now {ek}k=i,..., 2 u is the canonical basis in 
and for each initial condition we solve (34) in [0,Ty]. Let M be the transition matrix obtained 
by placing into adjacent columns the solutions Z{Ty) of (34) for all the initial conditions Z{{)) = 


{k = 1,2z/). Compute its eigenvalues and consider the largest in modulus: when such value is larger 
than 1 the longitudinal mode is torsionally unstable. 

In this section we implement the just described procedure and we describe the numerical computa¬ 
tions of the thresholds of torsional instability for each branch of periodic longitudinal modes. Several 
different indicators are used to measure the instability. The values of the parameters involved in the 


equations are taken as for the collapsed TNB, see Section 4.4 


Tablesummarises the thresholds of instability for the different branches of longitudinal modes: the 
second column is the energy threshold, the third column is the period of the longitudinal mode at the 
threshold, and the fourth is a measure of the amplitude of oscillation, that is, the gap between the 
maximum and the minimum of the longitudinal mode y\ 


A := maxy(x, t) 

x^t 


min^(x, t ). 

x.t 


Branch 

Energy (MJ) 

Period (s) 

A (m) 

1 

38. 

11.22 

5.8 

2 

51.8 

8.46 

10 

3 

15.5 

5.48 

6.5 

4 

53.7 

3.97 

7.8 

5 

74.1 

3.14 

6.9 

6 

56.6 

2.53 

4.5 

7 

91.4 

2.18 

5.2 

8 

95.8 

1.86 

4.6 

9 

87.1 

1.59 

3.8 

10 

82.1 

1.38 

3.3 


Table 2: Thresholds of instability on each branch of nonlinear modes. 


A careful look at the video m and the data in the Report [2] confirm that the oscillations prior 
to the TNB collapse were of the order of a few meters, as in our numerical results, see the values 
of A in Table The period at the threshold of instability should be compared with the period of 
the longitudinal mode on the same branch when the energy is close to 0, that is, when the nonlinear 
system (18) is close to linear, see Table 

Let us now introduce an effective tool to measure the torsional instability which allows the comparison 
of different longitudinal modes. The Floquet Theorem (see m Theorem 8.23]) states that, if T(t) 
is a fundamental matrix solution of (34), then there exists a matrix B and a T^-periodic matrix 
P such that T(t) = note that both B and P may be complex-valued. The eigenvalues 


Ai, X2iy of are called the characteristic multipliers of (34): if Ri,..., ¥ 2 ^ denote the corresponding 
normalized eigenvectors, then the solution Zj of ( [34| ) satisfying the initial condition Zj{0) = Vj (for 
some j — l,...,2z/) also satisfies Zj{Ty) = XjZj{0) and, in turn, Zj{mTy) — A^Zj(O) for any integer 
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m > 1. Whence represents the rate of growth of the amplitude of oscillation of Zj. This leads 

to the following definition: 


Definition 4. We call expansion rate of (34) the largest rate of growth for solutions of (34); 


STZ := max 

3 




For every branch 
the expansion rate 


of nonlinear longitudinal modes we compute the expansion rate of (34), that is, 
of the torsional component of the motion. Figure displays the expansion rate 



Figure 5: Expansion rate for modes 1 to 4 with respect to the energy. 


as a function of the energy for the branches of nonlinear longitudinal modes of order 1 to 4 whereas 
Figure [^displays the same for the modes of order 5 to 10. Note that the vertical axis (the ^^7^-axis) has 
different scales in the two figures. The reason is that £TZ is much larger for high nonlinear longitudinal 
modes. All the plots clearly show that there exists a threshold above which the torsional stability of 
the longitudinal mode is lost, that is, where the expansion rate exceeds 1. 

Table 1^ summarises the expansion rates on all the branches that we considered, for different values of 
the energy. Some values in Table are not available because the Newton algorithm would require too 
small increments of T in order to follow the branch of solutions. Note that the energy of a longitudinal 
mode is related to the amplitude of oscillation A and, for a given energy, higher modes have smaller 
amplitudes. Therefore, for the same amplitude of oscillation, the expansion rate of the 9th and 10th 
longitudinal mode is larger than that of the 5th, see also Table 
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Figure 6: Expansion rate for modes 5 to 10 with respect to the energy. 


Branch 

2 

4 

6 

8 

10 

12 

14 

1 

1 . 

1 . 

1.0662 

— 

— 

— 

— 

2 

1 . 

1.00365 

1.03904 

— 

— 

— 

— 

3 

1.00614 

1.02071 

1.02961 

1.08141 

1.20949 

— 

— 

4 

1 . 

1 . 

1.01287 

— 

— 

— 

— 

5 

1 . 

1 . 

1.00001 

1.01521 

1.50051 

— 

— 

6 

1 . 

1 . 

1 . 

1.09919 

1.16332 

— 

— 

7 

1 . 

1 . 

1 . 

1 . 

1.09852 

1.58567 

1.97158 

8 

1 . 

1 . 

1 . 

1 . 

1.00112 

1.66552 

— 

9 

1 . 

1 . 

1.01322 

1.01353 

1.24852 

1.76429 

2.12488 

10 

1 . 

1 . 

1 . 

1 . 

1.25447 

1.73715 

2.05263 


Table 3: Expansion rates versus energy (in MJ). 
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According to Eldridge [21 V-3], a witness on the day of the TNB collapse, the bridge appeared to be 
behaving in the eustomary manner and the motions were eonsiderably less than had oeeurred many 
times before. From [21 p.20] we also learn that in the months prior to the collapse one prineipal mode of 
oseillation prevailed and that the modes of oseillation frequently ehanged. On the day of the collapse, 
the torsional oscillations at the TNB started suddenly and, according to Farquharson [21 V-10], the 
motions, whieh a moment before had involved a number of waves (nine or ten) had shifted almost 
instantly to two. These observations show that the transfer of energy (from longitudinal to torsional) 
depends on the particular excited mode. 

Table [^ explains why torsional oscillations did not appear earlier at the TNB even in presence of 
larger longitudinal oscillations, as reported by Eldridge: there are longitudinal modes which have a 
very small expansion rate. Moreover, Table also explains why torsional oscillations appeared with 
smaller longitudinal oscillations: the 9th and 10th longitudinal mode appear more prone to generate 
torsional oscillations because they have large expansion rate STZ. 

4.4 Glossary and constants 

In this section we list the constants considered throughout the paper. 

• y — y{x,t) downwards vertical displacement of the center of the deck 

• 9 = 9{x,t) torsional angle of the deck 

• I half width of the deck (6m) 

• L length of the deck (853.44m) 

• M mass linear density of the deck (7198A;^/m) 

• E Young modulus of the deck and of the cables (210GPa) 

• I linear density of the moment of inertia of the cross section (0.15m^) 

• G shear modulus of the deck (81GPa) 

• K torsional constant of the deck (6.44 • 10“^m^) 

• m mass linear density of the cable (981A;^/m) 

• i^o horizontal component of the tension of the cables (58300A;A^) 

• H{x) tangential component of the tension of the cables 

• Lc length of the cable (868.62m) 

• So length of the longest hanger (72m) 

• A area of the section of the cable (0.1228m^) 

• Bridge at rest means bridge fully mounted, operational at equilibrium. 

• Bridge unloaded means that the deck is not mounted. 
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5 Conclusions 


For long time the aerodynamic effects have been considered the only cause of the TNB collapse. 
But, as quickly surveyed in the Introduction, no purely aerodynamic explanation is nowadays able to 
justify the origin of the torsional oscillation, which is the main culprit for the collapse. In this paper 
we suggest a new nonlinear model, which improves [5], and we focus our attention on the nonlinear 
structural behavior of suspension bridges. To this end, we follow a suggestion by Irvine [23] and 
we isolate the bridge from aerodynamic effects and dissipation. The new model and the resulting 
evolution PDFs which describe the dynamics of suspension bridges are presented in Section 

In Section we explain how to find periodic longitudinal oscillations as solutions of a nonlinear 
system of ODEs: contrary to linear systems, the period depends on the energy and on the amplitude 
of the oscillation. 

In Section [4T] we describe the theoretical framework governing the structural instability of suspension 
bridges. Our analysis is based on tools from the Floquet theory and a sufficient condition for the 
stability at low energy is obtained in Section [4^ Then, in Section 4.3, we put into the model the 


quantitative parameters of the collapsed TNB and we compute its thresholds of instability. The results 
are in line with the values observed on the day of the collapse and well fit the Report [2]. First, the 
computed critical amplitude of longitudinal oscillations is the same as in the video |12]. Second, our 
results explain why the TNB withstood larger longitudinal oscillations on low modes, but failed for 
smaller longitudinal oscillations on higher modes, see [21 pp.28-31]. Our results show that the origin 
of torsional instability is structural. 

It is clear that, in absence of wind or external loads, the deck of a bridge remains still. On the other 
hand, when the wind hits a bluff body (such as the deck of the TNB) the flow is modified and goes 
around the body. Behind the deck, or a “hidden part” of it, the flow creates vortices which are, in 
general, asymmetric. This asymmetry generates a forcing lift which starts the vertical oscillations of 
the deck. Up to some minor details, this explanation is shared by the whole community and it has 
been studied with great precision in wind tunnel tests, see e.g. [23 ESI iO]. The vortices induced 
by the wind increase the internal energy of the structure and generate wide longitudinal oscillations 
which look periodic in time. This is the point where our analysis starts, that is when the longitudinal 
oscillations of the structure reach a periodic motion which is maintained in amplitude by a somehow 
perfect equilibrium between the input of energy from the wind and internal dissipation. Throughout 
this paper we have seen that, in such situation, if the longitudinal oscillation is sufficiently large, 
then a structural instability appears: this is the onset of torsional oscillations. In order to obtain 
a more accurate description of the dynamics of suspension bridges, the next step should be to take 
into account also the behavior of the aerodynamic and dissipation effects after the appearance of the 
torsional oscillations. 
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